function plotrocks(filenumber)
    N=readdata(0,1);
%     S=readdata(0,2);
%     T=readdata(0,3);
%     multigridlevels=readdata(0,4);
%     bound=readdata(0,5);
time=readdata(0,6);
%     dtnom=readdata(0,7);
dx=readdata(0,8);
dz=readdata(0,9);
mx=readdata(0,10);
mz=readdata(0,11);
%     mcs=readdata(0,12);
%     mct=readdata(0,13);
%     mtype=readdata(0,14);
%     mTemp=readdata(0,15);
%     mstressxx=readdata(0,16);
%     mstressxz=readdata(0,17);
%     mstrainxx=readdata(0,18);
%     mstrainxz=readdata(0,19);
%     mfinitestrain=readdata(0,20);
%     vx=readdata(0,21);
%     vz=readdata(0,22);
%     P=readdata(0,23);
     Nsurfaces=readdata(0,24);
     surfaces=readdata(0,25);

mcs=zeros(1,N); mct=zeros(1,N);
xdivision=1; xplotcellsize=sum(dx)/xdivision;
x=0; for s=0:xdivision-1 ind=find(mx>=x & mx<x+xplotcellsize); mcs(ind)=s; x=x+xplotcellsize; end
zdivision=100; zplotcellsize=sum(dz)/zdivision;
z=0; for t=0:zdivision-1 ind=find(mz>=z & mz<z+zplotcellsize); mct(ind)=t; z=z+zplotcellsize; end
Ia=int32(find((-1).^mcs.*(-1).^mct>0));
Ib=int32(find((-1).^mcs.*(-1).^mct<=0));
clear('mcs','mct');

mtotalstrain=zeros(N,1); oldtime=time;

figure(2); clf; contstdhist=[]; conthist=[]; timehist=[];
figure(3); clf; hold on;
for count=filenumber
    
    N=readdata(count,1);
    S=readdata(count,2);
    T=readdata(count,3);
%     multigridlevels=readdata(count,4);
%     bound=readdata(count,5);
    time=readdata(count,6);
%     dtnom=readdata(count,7);
    dx=readdata(count,8);
    dz=readdata(count,9);
    mx=readdata(count,10);
    mz=readdata(count,11);
%     mcs=readdata(count,12);
%     mct=readdata(count,13);
    mtype=readdata(count,14);
    %mTemp=readdata(count,15);
    %mstressxx=readdata(count,16);
    %mstressxz=readdata(count,17);
    %mstrainxx=readdata(count,18);
    %mstrainxz=readdata(count,19);
    %mfinitestrain=readdata(count,20);
%     vx=readdata(count,21);
%     vz=readdata(count,22);
     %P=readdata(count,23);
     Nsurfaces=readdata(count,24);
     surfaces=readdata(count,25);

    dt=time-oldtime; oldtime=time;
    x=zeros(1,S); z=zeros(1,T); 
    for s=1:S-1 x(s+1)=x(s)+dx(s); end
    for t=1:T-1 z(t+1)=z(t)+dz(t); end
    figure(1); clf;
    
    avtemp=zeros(1,T);
    
    
%     [grid,xgrid,zgrid]=bilingrid(mtype,mx,mz,dx,dz,2,2);
%     [Sgrid,Tgrid]=size(grid); surface=zeros(1,Sgrid); moho=zeros(1,Sgrid);
%     surfaceind=zeros(1,Sgrid); mohoind=zeros(1,Sgrid); surface=zeros(1,Sgrid);
%     for s=1:Sgrid
%         depthslice=grid(s,:);
%         ind=find(depthslice<1);
%         [nothing,t]=max(zgrid(ind));
%         surface(s)=zgrid(ind(t));
%         surfaceind(s)=ind((t));
%         ind=find(depthslice<2.99);
%         [nothing,t]=max(zgrid(ind));
%         moho(s)=zgrid(ind(t));
%         mohoind(s)=ind(t);
%         ind=find(depthslice<0.5);
%         [nothing,indlow]=max(zgrid(ind)); indlow=ind(indlow);
%         ind=find(depthslice>=0.5);
%         [nothing,indhigh]=min(zgrid(ind)); indhigh=ind(indhigh);
%         a=(grid(s,indhigh)-grid(s,indlow))/(zgrid(indhigh)-zgrid(indlow));
%         b=grid(s,indlow)-a*zgrid(indlow);
%         surface(s)=(0.5-b)/a;       
%     end
subplot(3,1,2);
mTemp=readdata(count,15);
[grid,xgrid,zgrid]=bilingrid(mTemp,mx,mz,dx,dz,1,1); ind1=find(grid<1);
%[Zgrid,Xgrid] = meshgrid(zgrid,xgrid); pcolor(Xgrid,Zgrid,grid); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 500e3 0 300e3]); clear('Xgrid','Zgrid','mTemp');
contvals=[1100 1300 1400 1500]; cont=contourc(xgrid,zgrid,grid',contvals);
for t=1:T
    avtemp(t)=mean(grid(:,t));
end
%imagesc(xgrid,zgrid,grid'); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 500e3 0 300e3]); title([num2str(time/(60*60*24*365*1e6)) ' Ma ']); clear('mTemp');
[Zgrid,Xgrid] = meshgrid(zgrid,xgrid); grid(ind1)=nan; pcolor(Xgrid,Zgrid,grid); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 500e3 0 300e3]); title([num2str(time/(60*60*24*365*1e6)) ' Ma ']); clear('mTemp');
%print(gcf, '-djpeg100','-r600',['..\outfiles\temp' num2str(count)]); clf;

subplot(3,1,3); hold on;
sealevel=0.5*(surfaces(1,Nsurfaces)+surfaces(S,Nsurfaces));
surface=surfaces(:,1);
for sur=4:Nsurfaces
    if mod(sur,5)==0 plot(x,surfaces(:,sur)-sealevel,'k','markersize',0.5); end
end
plot(x,surfaces(:,sur)-sealevel,'g','markersize',2);
set(gca,'ydir','rev'); axis([0 500e3 -8000 15000]);
%print(gcf, '-djpeg100',['..\outfiles\surface' num2str(count)]); clf;

subplot(3,1,1);
mtype(Ib)=mtype(Ib)+0.5;
[grid,xgrid,zgrid]=bilingrid(mtype,mx,mz,dx,dz,2,2); ind0=find(grid<0.75);
[Zgrid,Xgrid] = meshgrid(zgrid,xgrid); grid(ind0)=nan; pcolor(Xgrid,Zgrid,grid); shading flat; set(gca,'ydir','rev'); axis equal; axis([0 500e3 0 300e3]); caxis([0 5.5]); colorbar; title([num2str(time/(60*60*24*365*1e6)) ' Ma ']); hold on; %clear('Xgrid','Zgrid','grid','mtype');


%imagesc(xgrid,zgrid,grid'); hold on;
%     for sur=1:Nsurfaces
%         plot(x,surfaces(:,sur),'k','markersize',2);
%     end
pos=1; allypoints=cell(1,length(contvals));
timehist(end+1)=time/(60*60*24*365*1e6);
while pos<length(cont)
    step=cont(2,pos);
    i=find(contvals==cont(1,pos));
    xcont=cont(1,pos+1:pos+step); ycont=cont(2,pos+1:pos+step);
    plot(xcont,ycont,'w');
    pos=pos+step+1;
    allypoints(i)={[cell2mat(allypoints(i)) ycont]};
end
for i=1:length(contvals)
    conthist(i,length(timehist))=mean(cell2mat(allypoints(i))); contstdhist(i,length(timehist))=std(cell2mat(allypoints(i)));
end
shading flat; set(gca,'ydir','rev'); axis equal; axis([0 500e3 0 300e3]); title([num2str(time/(60*60*24*365*1e6)) ' Ma ']); clear('grid','mtype');
%print(gcf, '-djpeg100','-r600',['..\outfiles\rock' num2str(count)]); clf;
print(gcf, '-djpeg100','-r600',['..\outfiles\rock_temp_surface' num2str(count)]);



figure(2); clf; hold on;
for t=1:length(timehist)
    for i=1:length(contvals)
        plot([timehist(t) timehist(t)],[conthist(i,t)-contstdhist(i,t) conthist(i,t)+contstdhist(i,t)]-4e4);
    end
end
for i=1:length(contvals)
    plot(timehist,conthist(i,:)-4e4,'r');
end

figure(3); hold on;
plot(avtemp,z); set(gca,'ydir','rev'); plot(adiabat(z,40e3),z,'r');
    
%     figure(3);
%     mstressxx=readdata(count,16);
%     mstressxz=readdata(count,17);
%     mstresstot=(mstressxx.^2+mstressxz.^2).^0.5; clear('mstressxx','mstressxz');
%     mstrainxx=readdata(count,18);
%     mstrainxz=readdata(count,19);
%     mstraintot=(mstrainxx.^2+mstrainxz.^2).^0.5; clear('mstrainxx','mstrainxz');
%     mvisc=mstresstot./(2*mstraintot); clear('mstresstot','mstraintot');
%     [grid,xgrid,zgrid]=bilingrid(mvisc,mx,mz,dx,dz,1,1); clear('mvisc');
%     imagesc(xgrid,zgrid,log10(grid')); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 500e3 0 300e3]); title([num2str(time/(60*60*24*365*1e6)) ' Ma ']); clear('mTemp');
%     print(gcf, '-djpeg100','-r600','visc'); clf;

    
%     [Sgrid,Tgrid]=size(grid); q=zeros(1,Sgrid); mohotemp=zeros(1,Sgrid);
%     for s=1:Sgrid %this loops assumes that the grid is the same used when finding the moho and surface
%         t=surfaceind(s);
%         q(s)=3*(grid(s,t+5)-grid(s,t+1))/(zgrid(5)-zgrid(1)); %assuming k=3
%         t=mohoind(s);
%         mohotemp(s)=grid(s,t);
%     end
%     clear('grid');
%     [ax,nothing,nothing]=plotyy(xgrid,q,xgrid,mohotemp); set(ax(1),'xlim',[0 1500000]); set(ax(2),'xlim',[0 1500000]); set(ax(1),'ylim',[3e-2 25e-2],'ytick',3e-2:1e-2:25e-2); set(ax(2),'ylim',[0 800],'ytick',0:50:800);
%     print(gcf, '-djpeg100',['..\outfiles\surfaceq-mohotemp' num2str(count)]); clf;

%     mstrainxx=readdata(count,18); mstrainxz=readdata(count,19);
%     mtotalstrain=mtotalstrain+dt*(mstrainxx.^2+mstrainxz.^2).^0.5;
%     [grid,xgrid,zgrid]=bilingrid(mtotalstrain,mx,mz,dx,dz,1,1);
%     %[Zgrid,Xgrid] = meshgrid(zgrid,xgrid); pcolor(Xgrid,Zgrid,grid); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 500e3 0 300e3]); clear('Xgrid','Zgrid','grid');
%     imagesc(xgrid,zgrid,grid'); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 500e3 0 300e3]); clear('grid');
%     print(gcf, '-djpeg100','-r600',['..\outfiles\totalstrain' num2str(count)]); clf;
%     
figure(1); clf;
mstrainxx=readdata(count,18);
mstrainxz=readdata(count,19);
[grid,xgrid,zgrid]=bilingrid(log10((mstrainxx.^2+mstrainxz.^2).^0.5),mx,mz,dx,dz,1,1); clear('mstrainxx','mstrainxz');
%[Zgrid,Xgrid] = meshgrid(zgrid,xgrid); pcolor(Xgrid,Zgrid,grid); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 500e3 0 300e3]); clear('Xgrid','Zgrid','grid');
[Zgrid,Xgrid] = meshgrid(zgrid,xgrid); grid(ind1)=nan; pcolor(Xgrid,Zgrid,grid); shading flat; set(gca,'ydir','rev'); axis equal; axis([0 500e3 0 300e3]); colorbar; title([num2str(time/(60*60*24*365*1e6)) ' Ma ']); hold on; clear('Xgrid','Zgrid','grid','mtype');
print(gcf, '-djpeg100','-r200',['..\outfiles\strain' num2str(count)]); clf;
% 
mfinitestrain=readdata(count,20);
[grid,xgrid,zgrid]=bilingrid(mfinitestrain,mx,mz,dx,dz,1,1); clear('mfinitestrain');
%[Zgrid,Xgrid] = meshgrid(zgrid,xgrid); pcolor(Xgrid,Zgrid,grid); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 500e3 0 300e3]); clear('Xgrid','Zgrid','grid');
[Zgrid,Xgrid] = meshgrid(zgrid,xgrid); grid(ind1)=nan; pcolor(Xgrid,Zgrid,grid); shading flat; set(gca,'ydir','rev'); axis equal; axis([0 500e3 0 300e3]); colorbar; title([num2str(time/(60*60*24*365*1e6)) ' Ma ']); hold on; clear('Xgrid','Zgrid','grid','mtype');
print(gcf, '-djpeg100','-r100',['..\outfiles\plasticstrain' num2str(count)]); clf;
    
%     mstressxx=readdata(count,16);    mstressxz=readdata(count,17);    mstrainxx=readdata(count,18);    mstrainxz=readdata(count,19);
%     P=readdata(count,23); mTemp=readdata(count,15); mtype=readdata(count,14);
%     zP=zeros(1,T)+0.5*dz(1); xP=zeros(1,S)+0.5*dx(1);
%     for t=1:T-2 zP(t+1)=zP(t)+0.5*(dz(t)+dz(t+1)); end
%     zP(T)=zP(T-1)+dz(T-1);
%     for s=1:S-2 xP(s+1)=xP(s)+0.5*(dx(s)+dx(s+1)); end
%     xP(S)=xP(S-1)+dx(S-1);
%     [ZP,XP] = meshgrid(zP,xP); P=P(1:end-1,1:end-1); XP=XP(1:end-1,1:end-1); ZP=ZP(1:end-1,1:end-1);
%     mP=interp2(XP',ZP',P',mx,mz);
%     stressinv=(mstressxx.^2+mstressxz.^2).^0.5; straininv=(mstrainxx.^2+mstrainxz.^2).^0.5;
%     
%     ind=find(mtype==3);
%     n=3.5; Ad=4.94e-21; E=360e3; F=power(2,(1-2*n)/n); F=0.5; R=8.314472; V=0;
%     %n=4; Ad=2.0000e-021; E=471e3; F=power(2,(1-2*n)/n); R=8.314472; V=15e-6;
%     expstress=2*F*Ad^(-1/n)*exp((E+V*mP(ind))./(n*R*(mTemp(ind)+273))).*straininv(ind).^(1/n);
%     loglog(expstress,stressinv(ind),'.'); hold on; loglog(expstress,expstress,'r.');
%     
%     ind=find(mtype==4);
%     expstress=2*F*Ad^(-1/n)*exp((E+V*mP(ind))./(n*R*(mTemp(ind)+273))).*straininv(ind).^(1/n);
%     loglog(expstress,stressinv(ind),'g.'); hold on; loglog(expstress,expstress,'r.');
%     
%     print(gcf, '-djpeg100',['..\outfiles\stresstest' num2str(count)]); clf;
    
    
% mtype=readdata(count,14); P=readdata(count,23); mstressxx=readdata(count,16); mstressxz=readdata(count,17); mTemp=readdata(count,15);
% zP=zeros(1,T)+0.5*dz(1); xP=zeros(1,S)+0.5*dx(1);
% for t=1:T-2 zP(t+1)=zP(t)+0.5*(dz(t)+dz(t+1)); end
% zP(T)=zP(T-1)+dz(T-1);
% for s=1:S-2 xP(s+1)=xP(s)+0.5*(dx(s)+dx(s+1)); end
% xP(S)=xP(S-1)+dx(S-1);
% [ZP,XP] = meshgrid(zP,xP); P=P(1:end-1,1:end-1); XP=XP(1:end-1,1:end-1); ZP=ZP(1:end-1,1:end-1);
% [gridP,xgrid,zgrid]=bilingrid(P(:),XP(:),ZP(:),dx,dz,1,1); clear ('XP','ZP','P');
% [grid,xgrid,zgrid]=bilingrid((mstressxx.^2+mstressxz.^2).^0.5,mx,mz,dx,dz,1,1); clear('mstressxx','mstressxz');
% [gridtype,xgrid,zgrid]=bilingrid(mtype,mx,mz,dx,dz,1,1); clear('mtype'); [Zgrid,Xgrid] = meshgrid(zgrid,xgrid);
% gridtype=round(gridtype); dens=zeros(size(gridtype));
% ind=find(gridtype<1); dens(ind)=1; ind=find(gridtype>=1); dens(ind)=2800; ind=find(gridtype>=2); dens(ind)=3000; ind=find(gridtype>=3); dens(ind)=3300;
% [gridTemp,xgrid,zgrid]=bilingrid(mTemp,mx,mz,dx,dz,1,1); clear('mTemp');
% dens=dens-3.25e-5.*gridTemp; clear('gridTemp');
% hydrostat=zeros(size(gridtype));
% for s=1:S
%     for t=2:T
%         hydrostat(s,t)=hydrostat(s,t-1)+9.82*0.5*(dens(s,t)+dens(s,t-1))*dz(t-1);
%     end
% end
% grid=grid-gridP+0*hydrostat; %sigmaxx=sigmaxxdev-(P-hydrastat)
% intgstress=zeros(1,S);
% for s=1:S
%     intgstress(s)=trapz(zgrid,grid(s,:));
% %     plot(zgrid,grid(s,:));
% %     plot(zgrid,gridP(s,:),'r');
% %     plot(zgrid,hydrostat(s,:),'g');
% %     plot(zgrid,gridP(s,:)-hydrostat(s,:),'c');
% %     plot(zgrid,grid(s,:)-gridP(s,:)+hydrostat(s,:),'k');
% end
% clear('grid','gridP','hydrostat','dens','gridtype');
% %clf; plot(intgstress);
% %print(gcf, '-djpeg100',['..\outfiles\stressdepth' num2str(count)]);
% avintgstress=trapz(xgrid,intgstress)/xgrid(end);
% figure(2); hold on; plot(time/(60*60*24*365*1e6),avintgstress,'r.'); figure(1);
% %pcolor(Xgrid,Zgrid,gridP-hydrostat); shading flat;
% set(gca,'ydir','rev'); axis equal; colorbar; axis([0 500e3 0 300e3]);

    
%     figure(1); plot(mx(I0(I0a)),mz(I0(I0a)),'y.','markersize',1); hold on; plot(mx(I0(I0b)),mz(I0(I0b)),'g.','markersize',1);
%     plot(mx(I1(I1a)),mz(I1(I1a)),'b.','markersize',1);
%     plot(mx(I1(I1b)),mz(I1(I1b)),'c.','markersize',1);
%     set(gca,'ydir','rev'); hold off; axis equal; title([num2str(time/(60*60*24*365*1e6)) ' Ma ']); axis([0 1000000 0 250000]);
%     print(gcf, '-djpeg100',['..\outfiles\a' num2str(count)]);
% 
%     [gridstrain,xgrid,zgrid]=gridmarkers(mfinitestrain,mx,mz,mcs,mct,x,z,S*2,T*2); [Xgrid,Zgrid] = meshgrid(xgrid,zgrid); figure(2); pcolor(Xgrid,Zgrid,gridstrain); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 1000000 0 250000]);
%     %[gridstrain,xgrid,zgrid]=gridmarkers1(mfinitestrain,mx,mz,mcs,mct,x,z); [Xgrid,Zgrid] = meshgrid(xgrid,zgrid); figure(2); pcolor(Xgrid,Zgrid,gridstrain); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 1000000 0 250000]);
%     print(gcf, '-djpeg100',['..\outfiles\plasticstrain' num2str(count)]); 
%     clear ('gridstrain','xgrid','zgrid'); clf;
%     I=find(mtype==0); mstrainxx(I)=0; mstrainxz(I)=0; %resetting air strain
%     [gridstrain,xgrid,zgrid]=gridmarkers((mstrainxx.^2+mstrainxz.^2).^0.5,mx,mz,mcs,mct,x,z,S*2,T*2); [Xgrid,Zgrid] = meshgrid(xgrid,zgrid); figure(5); pcolor(Xgrid,Zgrid,gridstrain); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 1000000 0 250000]);
%     %[gridstrain,xgrid,zgrid]=gridmarkers1((mstrainxx.^2+mstrainxz.^2).^0.5,mx,mz,mcs,mct,x,z); [Xgrid,Zgrid] = meshgrid(xgrid,zgrid); figure(5); pcolor(Xgrid,Zgrid,gridstrain); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 1000000 0 250000]);
%     print(gcf, '-djpeg100',['..\outfiles\' num2str(count)]);
%     clear ('gridstrain','xgrid','zgrid');
 
end
end